Get some US Census data

Download US Census data on natural resource dependent jobs (total # and median income) within counties and summarize across NEP boundaries:

nep_sf <- st_read(dsn = "./data-raw", layer = "NEP_Boundaries10162018", quiet = TRUE)
states <- c("AL", "AZ", "AR", "CA", "CO", "CT", "DE", "FL", "GA",
            "ID", "IL", "IN", "IA", "KS", "KY", "LA", "ME", "MD",
            "MA", "MI", "MN", "MS", "MO", "MT", "NE", "NV", "NH", "NJ",
            "NM", "NY", "NC", "ND", "OH", "OK", "OR", "PA", "RI", "SC",
            "SD", "TN", "TX", "UT", "VT", "VA", "WA", "WV", "WI", "WY",
            "PR")

metrics <- load_variables(2017, "acs5", cache = TRUE) %>% 
       filter(grepl("Fishing|fishing", label))

totaljobs_sf <- get_acs(geography = "county", variables = c(Ag_For_Fish_Hunt_Min = "C24050_002"), 
                         state = states, geometry = TRUE)

nr_medinc_sf <- get_acs(geography = "county", variables = c(Med_Income = "B24031_003"), 
                         state = states, geometry = TRUE)

ustotaljobs_sf <- st_transform(totaljobs_sf, st_crs(nep_sf))
usnrmedinc_sf <- st_transform(nr_medinc_sf, st_crs(nep_sf))

nep_jobs_intersects <- st_intersects(nep_sf, ustotaljobs_sf)
nep_medinc_intersects <- st_intersects(nep_sf, nr_medinc_sf)

nep_sel_sf <- ustotaljobs_sf[unlist(nep_jobs_intersects),]
nep_sel2_sf <- nr_medinc_sf[unlist(nep_medinc_intersects),]

nep_jobs <- st_join(nep_sf, nep_sel_sf, join = st_intersects) %>%
                   sf::st_buffer(dist = 0)
nep_nrmedinc <- st_join(nep_sf, nep_sel2_sf, join = st_intersects) %>%
                   sf::st_buffer(dist = 0)

nep_jobs_sum <- nep_jobs %>% 
               select(NEP_NAME, estimate) %>%
               group_by(NEP_NAME) %>% 
               summarise(jobs = sum(estimate))

nep_medinc_sum <- nep_nrmedinc %>% 
               select(NEP_NAME, estimate) %>%
               group_by(NEP_NAME) %>% 
               summarise(medinc = median(estimate))

US Census Natural Resource Dependent Jobs by NEP Project Areas

Overview Map

Plot US Census natural resource dependent jobs data within the NEP boundaries:

pal <- colorNumeric(palette = "viridis", domain = nep_jobs_sum$jobs, n = 10)

nep_jobs_sum %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Positron") %>%
    addPolygons(popup = ~ paste(NEP_NAME,"</br>","Jobs = ",jobs),
                stroke = FALSE,
                smoothFactor = 0,
                fillOpacity = 0.7,
                color = ~ pal(jobs)) %>%
    addLegend("bottomright", 
              pal = pal, 
              values = ~ jobs,
              title = "Natural Resource Dependent Jobs</br>(Ag.,Forest.,Fishing,Hunting,Mining)",
              opacity = 1)

Data Table

Underlying natural resource dependent jobs estimated within NEP Project Areas:

nep_jobs2 <- nep_jobs_sum  %>% 
             select(NEP_NAME, jobs) %>% 
             st_set_geometry(NULL) %>% 
             rename("NEP" = NEP_NAME, "Total Jobs(2015)" = jobs)

nep_jobs2

US Census Natural Resource Dependent Jobs, Median Incomes by NEP Project Areas

Overview Map

Plot US Census natural resource dependent jobs, median income data within the NEP boundaries:

pal2 <- colorNumeric(palette = "viridis", domain = nep_medinc_sum$medinc, n = 10)

nep_medinc_sum %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Positron") %>%
    addPolygons(popup = ~ paste(NEP_NAME,"</br>","Median Income = ",medinc),
                stroke = FALSE,
                smoothFactor = 0,
                fillOpacity = 0.7,
                color = ~ pal2(medinc)) %>%
    addLegend("bottomright", 
              pal = pal2, 
              values = ~ medinc,
              title = "Natural Resource Dependent Jobs</br>(Median Annual Income)",
              opacity = 1)

Data Table

Underlying natural resource dependent jobs, median incomes estimated within NEP Project Areas:

nep_medinc2 <- nep_medinc_sum  %>% 
             select(NEP_NAME, medinc) %>% 
             st_set_geometry(NULL) %>% 
             rename("NEP" = NEP_NAME, "Median Income ($, 2015)" = medinc)

nep_medinc2

US BEA Natural Resource Dependent Jobs, Total Personal Income by NEP Project Areas

Import US BEA data on natural resource dependent jobs (total county income) and summarize across NEP boundaries:

userSpecList <- list('UserID' = '88964E07-8671-4DFB-BF65-4E680DEDA641',
                     'Method' = 'GetData',
                     'datasetname' = 'Regional',
                     'TableName' = 'CAINC5N',
                     'LineCode' = '100',
                     'GeoFIPS' = 'COUNTY',
                     'Year' = '2017')

nr_income <- beaGet(userSpecList, asTable = TRUE) %>% 
                    filter(str_detect(GeoFips, "....0")==FALSE) %>%    
                    mutate(GEOID = as.character(GeoFips))
nr_income$GEOID <- gsub(" ", "", nr_income$GEOID, fixed = TRUE)

nr_inc_sf <- left_join(nep_sel_sf, nr_income, by = c('GEOID' = 'GEOID'))

nep_totinc <- st_join(nep_sf, nr_inc_sf, join = st_intersects) %>%
                   sf::st_buffer(dist = 0)

nep_inc_sum <- nep_totinc %>% 
               select(NEP_NAME, DataValue_2017) %>%
               group_by(NEP_NAME) %>% 
               summarise(totinc = sum(DataValue_2017, na.rm = TRUE))
nep_inc_sum$totinc <- ifelse(nep_inc_sum$totinc==0,NA,nep_inc_sum$totinc*1000)

Overview Map

Plot US BEA natural resource dependent jobs total income estimates within the NEP boundaries (where available):

pal3 <- colorNumeric(palette = "viridis", domain = nep_inc_sum$totinc, n = 10)

nep_inc_sum %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Positron") %>%
    addPolygons(popup = ~ paste(NEP_NAME,"</br>","Total Income (2017) = ",totinc),
                stroke = FALSE,
                smoothFactor = 0,
                fillOpacity = 0.7,
                color = ~ pal3(totinc)) %>%
    addLegend("bottomright", 
              pal = pal3, 
              values = ~ totinc,
              title = "Natural Resource Dependent Jobs</br>(Total Annual Income, if reported)",
              opacity = 1)

Data Table

Underlying natural resource dependent jobs, total personal income estimated within NEP Project Areas:

nep_inc2 <- nep_inc_sum  %>% 
             select(NEP_NAME, totinc) %>% 
             st_set_geometry(NULL) %>% 
             rename("NEP" = NEP_NAME, "Total Personal Income ($, 2017)" = totinc)

nep_inc2